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Abstract 

Baryons as relativistic bound states in 3-quark correlations are described by an ef- 
fective Bethe-Salpeter equation when irreducible 3-quark interactions are neglected 
and separable 2-quark correlations are assumed. We present an efficient numeri- 
cal method to calculate the nucleon mass and its covariant wave function in this 
quantum field theoretic quark-diquark model with quark-exchange interaction. Ex- 
panding the components of the spinorial wave function in terms of Chebyshev poly- 
nomials, the four-dimensional integral equations are in a first step reduced to a 
coupled set of one-dimensional ones. This set of linear and homogeneous equations 
defines a generalised eigenvalue problem. Representing the eigenvector correspond- 
ing to the largest eigenvalue, the Chebyshev moments are then obtained by iteration. 
The nucleon mass is implicitly determined by the eigenvalue, and its covariant wave 
function is reconstructed from the moments within the Chebyshev approximation. 
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PROGRAM SUMMARY 



Title of program: BSE 

Computers: Workstation DEC Alpha, LINUX PCs (AMD K7) 

Operating system under which the program has been tested: UNIX, Linux 

Programming language used: Fortran 90 

Memory required to execute with typical data: 10 MB 

No. of bits in a word: 32 

Peripherals used: standard output, disk 

No. of lines in distributed program, including test data, etc.: 

Keywords: linear integral equations; Bethe-Salpeter equation; nucleon model; 
diquarks. 

Nature of physical problem: Diquarks are introduced as separable correlations 
in the two-quark Green's function to facilitate the description of baryons as 
relativistic three-quark bound states. These states then emerge as solutions 
of Bethe-Salpeter equations for quarks and diquarks that interact via quark 
exchange. 

Method of solution: Chebyshev polynomials are used for an expansion of the 
Bethe-Salpeter vertex and wave functions and an (approximative) represen- 
tation of the interaction kernel and propagators. The resulting set of coupled 
one-dimensional integral equations for the corresponding moments of vertex 
and wave function is converted to a matrix equation by employing Gaussian 
quadrature grid points. This equation is then solved iteratively for a certain 
test mass of the bound state. A bisection method in the bound state mass 
variable is used to determine the correct mass corresponding to the chosen 
quark and diquark parameters. 

Restrictions on the complexity of the problem: The separable quark-quark cor- 
relations are restricted to the leading covariants in the scalar and axialvector 
(i.e. positive parity) channels. 

Typical running time: Depending on the parameters and the desired accuracy 
between one and ten minutes on 1.2 GHz AMD K7. 
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LONG WRITE-UP 



1 The physical problem 

A Poincare- invariant, quantum field theory based model of the nucleon is re- 
quired to describe the high-precision data of nucleon properties in the several- 
GeV energy range. To this end, we here demonstrate how to calculate (using 
certain approximations detailed below) the matrix elements of three quark 
operators q between the vacuum \Q) and the nucleon bound state \Pn), the 
asymptotic momentum state with nucleon quantum numbers. This matrix el- 
ement ip ~ (n\T(qqq)\Pfir) is the covariant Bethe-Salpeter wave function of 
the nucleon which can be applied as input in covariant calculations of many 
observables such as the various form factors of the nucleon. In the following sec- 
tions we describe an approximation scheme based on the Faddeev equations 
(see, e.g., ref. [1]) which makes use of the concept of diquarks as separable 
quark-quark correlations [2-8]. In this scheme, the covariant diquark-quark 
model, see refs. [9-16], we obtain a tractable set of equations for the com- 
ponents of the covariant nucleon wave function which can be solved without 
non-relativistic reductions. 

1.1 The separable approximation to the relativistic 3- quark problem 

To summarize the relation between the covariant quark-diquark model and 
the general, relativistic three-quark problem briefly in the following, we start 
from the six-quark Green function, 

G(x i ,y i ) = (0\Tf[q(x i )q(y i )\0). (1) 

i=i 

The variables Xi and yi not only represent the space-time coordinates of the 
quark fields, but also include their discrete labels for color, spin, and flavor. 
The six-point function (1) satisfies the Dyson equation 

G = G + G o K o G . (2) 

In this equation the disconnected six-point function G describes the free 
propagation of three dressed quarks, and the three-quark scattering kernel K 
contains all two- and three-particle irreducible diagrams. The symbol "o" in 
eq. (2) denotes summation/integration over all independent internal coordi- 
nates and labels which defines a composition law for linear maps on a suitable 
function space. Unless explicitly stated otherwise we will henceforth work in 
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momentum space with Euclidean metric. It is thus not necessary to introduce 
different symbols for momentum and coordinate space objects. 



The nucleon as a three-particle bound state with mass M manifests itself as 
a pole in the six-point function at — P 2 = M 2 where P = p 1 + p 2 + Pz is the 
total four-momentum of the three-quark system. One can thus parameterize 
the six-point function in the vicinity of the pole as 

nn \ ^(k 1 ,k 2 ,k 3 ) $(pi,p 2 ,p 3 ) , Q x 
G(h, Pi) P 2 + M 2 > ( 3 ) 

where denotes the bound state wave-function. Substituting this parameteri- 
zation into the Dyson equation (2) and identifying the residuum contributions 
on both sides, one obtains the homogeneous bound state equation, 

i> = G o K o ^ G' 1 o ijj = . (4) 

Despite the simple appearance of the bound state equation in this symbolic 
form, its direct solution is not feasible. Even in principle, this could only be 
attempted once a model interaction is specified to fix the contributions to the 
kernel K which can in general never be fully determined from first principles. 

The approximation scheme we present below essentially consists of two steps. 
Thereby, some of the contributions to the kernel, namely those due to irre- 
ducible three-particle interactions, are neglected while an explicit reference to 
the others is avoided by shifting the details of the two-particle interactions 
into the parametrisations of the relevant diquark correlations. 

The first step leads to the corresponding relativistic Faddeev equations. Ne- 
glecting all contributions from irreducible three-quark interactions, the kernel 
K can be written as a sum of three two-quark interaction kernels, 

K = K x + K 2 + K 3 . (5) 

We adopt the notation that the subscript of _£Q refers to the spectator quark 
qi. The Ki = k qq <E> S^ 1 are defined in 3-quark space and thus contain an 
inverse quark propagator S^ 1 acting on the subspace of the spectator quark. 
The quark pair interacting via k qq in Ki then is (qj,qk) with the three labels 
k) being a cyclic permutation of (1, 2, 3). These 2-quark kernels determine 
the interactions for disconnected scattering amplitudes Tj = t qq <g> S*^ 1 which 
describe the scattering between quark j and k (with amplitude t qq ) and which 
are summed in the Dyson series for this 2-quark subspace as follows: 

T i = K i + K i oG oT i . (6) 

Introducing the Faddeev components of the 3-quark wave function ip, 

i>i = G oKioi>, (7) 
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one then has — i>i by virtue of eqs. (4,5). Furthermore, eq. (6) can now be 
used to simplify eq. (4) considerably yielding the coupled Faddeev equations, 

ipi = G o T t o (tpj + fa) , (8) 

for Faddeev components with pairwise distinct k G {1, 2, 3}. 

The assumption central to the second step, and the defining one of the quark- 
diquark model, is the separability of the two-quark scattering amplitude t qq . 
This amplitude is thereby typically parametrised by scalar and axialvector di- 
quark correlations to account for its most important separable contributions 
at sufficiently small (absolute values of the complex Euclidean) diquark mo- 
menta, 



t qQ (k j ,k k ; Pj , Pk )=x 5 ((kj - h)/2) D 55 (kj + k k ) f((p 3 - Pk )/2) + 

■xTUJzj - h)/2) D^(kj + k k ) x v {{Pj ~ Pk)/2) ■ (9) 

Here we employ simple particle-pole contributions to represent these diquark 
correlations. The diquark propagators are then given by the free scalar prop- 
agator and the Proca propagator for a spin-1 particle, respectively, 



D 55 (p) = 
D^(p) = 



1 

p 2 + m 2 sc 
1 

p 2 + m 2 ax 



m 2 



(10) 
(11) 



The diquark-quark vertices x i n e Q- (9) need to be antisymmetric with re- 
spect to the interchange of their quark indices. Both diquarks belong to the 
antisymmetric color antitriplet. The scalar diquark vertex is an antisymmetric 
isospin singlet and is also antisymmetric in Dirac space while the axialvector 
diquark vertex is symmetric in both Dirac and isospin indices (triplet). We do 
not give their color and flavor structures explicitly here, and we restrict the 
Dirac structures of both x 1° the dominant ones which are given by 



X%(p)=9s(l 5 C) jk V s (p), (12) 
xSk(p)=9a(-fC)ikV a (p), (13) 

respectively, where g s and g a represent coupling constants which are deter- 
mined by the diquark normalizations, see below. The scalar functions V s>a (p) 
parametrise the extensions of the two diquark-quark amplitudes in momentum 
space. Phenomenologically, a dipole form for their dependence on the quark 
relative momenta has proven successful, and we will assume equal widths for 
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both diquarks here, 

= = (377) < 14 > 

to reveal the general structures in a reasonably simple form in the following. 

The diquark widths determine their couplings g s and g a by normalization con- 
ditions derived from the Bethe-Salpeter equation (6) for the 2-quark scattering 
amplitude, 

tqq = k qq + k qq O (S <g) S) O t qq . (15) 

Inserting the ansatz (9) with the diquark pole contributions (10,11) into in the 
inhomogeneous equation (15) implicitly fixes the normalizations adopted for 
the diquark amplitudes. These implicit conditions are most easily derived from 
the derivative of eq. (15) w.r.t. the total diquark momentum P d = kj + k k , 
for details, see e.g. refs. [9,14]. Under the mild additional assumption that 
the 2-quark interaction kernel k qq be independent of P^, the normalization 
conditions for the scalar and axialvector diquark amplitudes are, 

4mL = t o (p d ■ ^-(S ® S)) o x 5 , (16) 

12mL = x"o(p < ,~(5®5))o^ (^ + ^?) , (17) 

respectively. The additional factor of 3 in the last line arises from the sum 
over the polarisation states of the spin-1 particle. 

1.2 The quark-diquark Bethe-Salpeter equation 

For identical particles, the three Faddeev components (7) are obtained from 
a unique amplitude 0^ by cyclic permutations of the indices and arguments 
labelling the individual quarks, and the 3-quark amplitude is obtained by their 
sum or, equivalently, their Faddeev wave function by ip = GooJ2 cy ciic 'Pijk- With 
the separable form for the two-quark scattering amplitude t qq in eq. (9), the 
following ansatz is employed in the Faddeev equations (8): 

<Pijk(Pi,Pj,Pk)= (18) 
X a jk (q)D 3b ((l-rj)P-p) (<S> b (p,P)u(P)). : (a, b = 1 . . . 5) . 

Herein, u(P) is a positive-energy spinor for the nucleon state with momentum 
P = Pi +pj +Pk, P = (l—r))pi — r]{pj +Pk) is the relative momentum between 
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p = X\P+p 




P d =(l-r\)P-p -p-k+(l-2T])P = q 

Fig. 1. The Bethe-Salpeter equation for the vertex function <J>. 

quark and diquark, and q = (pj — Pk)/2 is the relative momentum of the two 
quarks within the diquarks. We have introduced a momentum partitioning 
parameter 77 e [0, 1] which distributes the total momentum P between quark 
and diquark. An analogous parameter could also be introduced in a more gen- 
eral definiton of the relative momentum q within the diquarks as discussed 
in [14]. Translational invariance entails that observables, i.e. matrix elements 
like the form factors etc., should not depend on such momentum partition- 
ing parameters in the covariant bound-state wave functions. Furthermore, in 
eq. (18), $ a is composed of Dirac matrices such that $ 5 ^-u is the most general 
effective [vector] spinor compatible with nucleon quantum numbers for the 
Faddeev components. We call $ a the nucleon Bethe-Salpeter vertex function 
with quark and diquark. 



The Faddeev equations (8) with the ansatz (18) now reduce to a system of ef- 
fective 2-particle Bethe-Salpeter equations for a quark-diquark state bound by 
repeated quark-exchanges. Introducing an eigenvalue A, for the lowest bound 
state in the given channel, these equations have to be solved for the largest 
value of A = A(P 2 ) at a given P 2 , 



A(P 2 ) $ a (p, P)=J-0r 4 K ab (p, k, P) ^ b (k, P) , with (19) 
ty a (p,P) = S(r]P + p)D ab ((l-r])P -p) $ b (p,P) , and (20) 



K a >,A;,P) = -- 



-^X 5 (pi) S T ( q ) x 5 (p 2 ) # x"(pi) sT M X 5 (P2 



\ 4 xV) s T ( q) x a ( P 2) ^xV) s T ( q) x a ( P 2) 



1 



, (21) 



with the additional constraint that \(P 2 ) = l/g 2 at P 2 = M 2 , thereby implic- 
itly determining the nucleon mass. We introduced the Bethe-Salpeter wave 
function ^ = Go o $ by attaching quark and diquark legs to the vertex func- 
tion. With (20) in eq. (19) this determines a generalized eigenvalue problem 
in which the vertex function $ = K o Go o $ represents an eigenvector of the 
kernel K with respect to a metric given by the disconnected quark-diquark 
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propagator, 

G ab = S(rjP + p) D 3b ((l-rj)P-p) (2tt) 4 5\p - k) . (22) 

The numerical factors in the elements of the quark-exchange kernel K ab arise 
from the structure of the color and flavor couplings in the Faddeev equations. 
Momentum conservation fixes the momentum of the exchanged quark, 

q = -p-k + (l-2 V )P , (23) 

and the relative momenta of the outgoing and the incoming quark-pair in 
x{pi) and x{P2)-> respectively, are (see Fig. 1), 



Pl =p + k/2-(l-3r))P/2, (24) 
p 2 = -k-p/2 + (1 -3r])P/2 . (25) 

Since neither the eigenvalue A = l/g 2 s nor the mass M n of the bound nucleon 
are allowed to depend on the momentum partitioning parameter 77, for every 
solution ^/(p,P;r]i) of the Bethe-Salpeter equation (19) there should exist a 
whole family of solutions of the form ^/(p + (772 — Vi)P> P\ ^2)- 



1.3 Decomposition of wave and vertex function 



Having defined the effective spinor of the spectator quark to be $ a w in eq. (18) 
guarantees that the 3-quark wave function ip describes a spin-1/2 object. 
Furthermore, $ a can most generally be chosen to be an eigenfunction of 
A + = 1 + If I (iM), the positive energy projector, since A~u = (1 — A + )u = 0. 
Requiring positive parity for the Faddeev amplitudes leads to the additional 
condition 



7 4 $ 5 (p, P) 7 4 



$ 5 G»,p) 



\^(p,p)hy \-* v (p,p) 



(26) 



with the parity operation on a 4-vector defined by q = (—q, q A ) T . 



The positive-energy constraint and condition (26) greatly restrict the number 
of independent components in the vertex function. The scalar correlations 
$ 5 are described by two components and the axialvector correlations by six 
components, 



U\p,p)\ 




\^(p,p)J 





I 



ES t (p 2 ,p-P)S l (p,P) 



\ 



i=l 



E A{p 2 ,P- P) i 5 A?( P ,P) 



(27) 



8 



The scalar functions Si and Ai depend on the two independent Lorentz in- 
variants (p 2 and p ■ P) that can be formed from the relative momentum p 
and the on-shell nucleon momentum P. Normalized 4-vectors are introduced 
as p = p/\p\, and for the special case of the complex (i.e., timelike) nucleon 
momentum we adopt the convention that P = P/(iM), see below. 

The Dirac components describing the scalar correlations may be built out of 
A + and p"A + and the Dirac part of the axialvector correlations can be con- 
structed using the matrices P M A + , P^p ! A + , 7 M A + , 7 M ^A + , p^A + and p^p'A + . 
Herein, for later convenience, we select a decomposition such that the scalar 
functions Si and Ai decouple within the Dirac and Lorentz components of (27) 
in the nucleon rest-frame, 



(0, iMf 



(28) 



and with additionally choosing the spatial part of the relative momentum to 
point in the z-direction (see sec. 2.1), 

p = (0, 0, \p\ simp, \p\ cos^) T . (29) 

A decomposition of the covariants achieving this turns out to be given by 



51 = A+ 

5 2 = -iJ T A+ 



At 



At 




At 


= P»A+ 


At 


= ptPrA + 


At 


= j*£A+ 


At 


= ^ T A+-At 


,At 


= ^ T i> T A+-A 



(30) 



Here and in the following, the subscript t denotes contraction with 5^ v — P^P V 
to project 4-vectors onto the subspace transverse to the nucleon momentum, 
e.g., px = p—P(p-P). While the decomposition (30) is well suited for numerical 
processing, there exists an alternative decomposition of the Faddeev amplitude 
into components which for a nucleon at rest lead to spin and orbital angular- 
momentum eigenstates [9]. The uniquely defined amplitudes in this partial 
wave decomposition are certain linear combinations of the functions Si and 
Ai as defined by eqs. (27,30) above. 

Since the constraints on the vertex function apply in the same way also to 
the quark-diquark Bethe-Salpeter wave function \1> defined in eq. (20), its 
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analogous decomposition reads, 



W(p,P) 



ts i (p\p-P)S t {p,P) S 

i=l 



\i=l I 



(31) 



The coefficient functions Si, Ai herein can be expressed in terms of the Si, Ai 
that occur in the decomposition (27) of the vertex function $. The explicit 
expressions follow readily from the relation \I> = G o $, c.f., eqs. (20,22). 

With the formalism being fully Poincare covariant we are free to choose the 
Lorentz frame that best suits our purpose of solving the Bethe-Salpeter equa- 
tions (19-21). As a particularly convenient choice we adopted the nucleon rest 
frame for the calculations presented herein, with the total and relative mo- 
menta, P and p, as given in eqs. (28) and (29), respectively. 

The Bethe-Salpeter equations (19-21) are rewritten in terms of equations for 
the scalar functions S±, . . . , A 6 . Their numerical solution described in the fol- 
lowing section can be divided into two major steps: 

• Chebychev expansion of these scalar functions to account for their p ■ P 
dependences and Chebyshev polynomial approximation of the product of 
the free quark and diquark propagators in the variable p ■ P and of the 
quark exchange kernel in the two variables p ■ P and k ■ P. 

• Solution by numerical iteration of the coupled system of one- dimensional 
integral equations resulting for the Chebychev moments of the scalar func- 
tions. 



Note that the expansion in terms of Chebychev polynomials employed in step 
two, though similar in nature, is not quite the same as a hyperspherical ex- 
pansion. The Chebychev expansion turns out to be extremely efficient here. 
Since typically only a few Chebychev polynomials are needed, see below, the 
numerical effort is greatly reduced as compared to any attempt at calculating 
the amplitude and the wave function directly on a two-dimensional grid. 

The underlying reason for this efficiency relates to the symmetries of Bethe- 
Salpeter equations in certain limits. For a more detailed discussion of the 
analogous method within the Wick-Cutkosky model see, e.g., section one of 
ref. [17] and the references therein. The numerical solution of the ladder- 
approximate Bethe-Salpeter equation in the massless as well as the massive 
Wick-Cutkosky model has been studied extensively in the literature, for a 
review see [18] and for a detailed numerical investigation see ref. [19]. This 
model describes the interactions of massive scalar particles by a likewise scalar 
but possibly massless particle exchange. In the case of a massless exchange- 
particle, the Bethe-Salpeter equation exhibits an 0(4) symmetry [20]. By 
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expanding the wave functions into hyperspherical harmonics one finds that 
the zero orbital angular momentum states of the Wick-Cutkosky model are 
essentially given by Chebyshev polynomials of the second kind U n (z). These 
are related to the Chebyshev polynomials of the first kind T n (z), used in our 
expansions below, by [21] 

T n (z) = U n (z)-zU n - 1 (z), (n>0). (32) 

This thus exemplifies the close relation between the expansions of wave func- 
tions in terms of the T n 's and those into hyperspherical harmonics. 

The quark-diquark Bethe-Salpeter equations (19-21), as derived from the rel- 
ativists Faddeev equations under the separable assumption, result to be of a 
ladder-exchange type also. Here, the mixing of different angular momentum 
components leads to additional structures, however. This is due to the spinor 
nature of the Bethe-Salpeter wave function. Much the same as for ordinary 
Dirac spinors, its lower components have different angular momentum than 
the upper ones in the relativistic treatment of a spin-| bound state. In or- 
der to obtain a closed system of equations, different partial waves have to 
be taken into account. The expansion into spinor hyperspherical harmonics is 
nevertheless possible and it convergences rapidly due to an approximate 0(4) 
symmetry for fermion-boson Bethe-Salpeter equations of the kind under con- 
sideration [14,15]. It is therefore understandable also for physical reasons that 
the expansion in terms of Chebychev polynomials provides such an efficient 
numerical technique. 

To conclude this section, we discuss the kinematics which can be employed to 
test the procedure and its limitations. In principle, the eigenvalue A = 
should be independent of the Mandelstam parameter rj (see the discussion 
after eq. (21)). Due to singularities in the propagators, however, without con- 
siderable modifications to the numerical procedure, it is not possible to use the 
whole interval [0,1] for such a test [9,14,15]. Rather, the range allowed for the 
momentum partitioning parameter 77 in the nucleon Bethe-Salpeter equation 
is restricted to 



rj G [1 - m sc /M n , m q /M n ] , (33) 

if m ax > m sc > m q is assumed. Singularities in the exchange quark propagator 
and the diquark vertices employing n-pole scalar functions V, see eq. (14), lead 
to the additional bounds 



77 G - m q /M n ), |(1 + m q /M n )\ , (34) 
77 G [1(1 - 2X n /M n ), |(1 + 2X n /M n ] . (35) 

Note that these bounds on the value of 77 arise in the practical calculations 
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when performed as outlined above. Though an arbitrary rj in [0,1] could be 
chosen in principle, additional residue terms have to be included beyond these 
bounds in the Euclidean formulation to account for its proper connection with 
the underlying Bethe-Salpeter equation in Minkowski space. 



Well within the above bounds on the momentum routing, however, pronounced 
plateaus can be observed in order to verify the 77-independence of physical 
observables [9,14]. The extension of these plateaus increases with increasing 
orders in the Chebyshev expansions. 

As the values of 77 approach the boundaries given in eqs. (33-35), however, 
strong variations occur due to the vicinity of the poles in the propagators with 
the effect that the convergence of their Chebyshev expansion, c.f., eq. (47) be- 
low, is slowed down considerably. In order to achieve some sufficient accuracy 
on the vertex function, larger and larger numbers of Chebyshev moments n max 
are then needed to approximate the components of the wave function if). The 
slower convergence of the wave function expansion under these circumstances 
is also observed in the numerical solutions [9,14,15]. 



2 Numerical method 



2. 1 Projection on scalar functions in nucleon rest frame 



In the nucleon rest frame (28,29) we express the Euclidean components of p^ 
and k^ (the latter being the integration variable) in hyperspherical coordi- 
nates: 



k 2 
k 3 



\k\ 



^ sin ip' sin 9' sin 0'^ 
sin ip' sin 9' cos <j)' 
sin ip' cos 9' 

y cos ip' j 



\P 4 J 



\P\ 



< ^ 


sin-0 
\^cos ip ) 



(36) 
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To simplify the notation, we set z = cosip, z' = cosip'. In the chosen frame 
the vertex function decomposition on the l.h.s of eq. (19) simplifies to 



$3 



{p\z) = 



( 



Si(p 2 ,z) 







\a 3 VT^S 2 (p 2 ,z) 
^sVT^A^p^z) 
v A 2 (p 2 ,z) 
' ia 3 A 3 (p 2 ,z) O" 
K iVT^A 4 {p 2 ,z) Oj 
' ia 2 A 5 {p 2 ,z) 
K -a, VT^A 6 (p 2 ,z) 
' %a x A^p 2 ,z) o\ 
v \a 2 ^l- z 2 A % {p 2 ,z) OJ ) 



(37) 



By virtue of the positive-energy condition and our choice for the covariants in 
eqs. (30) the vertex (and likewise the wave function) consists of 2x2-blocks 
in Dirac space which define upper and lower components for the scalar and 
Lorentz components of $, respectively. The unknown scalar functions do not 
couple within these blocks. Since we have chosen the most general decompo- 
sition, the equations for the Lorentz components $ 2 and are degenerate. 

For convenience we will use generic functions Yi {% — 1, . . . , 8) related to the 
functions Si and A^ 



5*1,2 —> Yi' 



A,. 



Y 



3.. .8 • 



(3* 



The functions Y { substituting the functions Si, A\ appearing in the wave func- 
tion (31) are defined analogously. 

Projection onto the scalar functions Yi and Yi is now done by inspecting the 
upper and lower components in the equations (20,19) for each Lorentz com- 
ponent. As explained before, the Yi and Yi neatly decouple in upper and lower 
components. We arrive at: 



Y t (p 2 ,z) = (g^(p 2 ,z)Y 3 (p 2 ,z), (39) 
Y l {p 2 ,z)=j-^{H'f{p 2 ,k 2 ,y',z,z')Y l {k 2 ,z'). (40) 

As a result, we obtain the propagator matrix (<?o) u an d the modified quark 
exchange kernel matrix (H')^ . The latter depends on the possible scalar prod- 
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ucts between the vectors k,p,P. These scalar products can be expressed as 
k 2 , p 2 , z' = k ■ P, z = p ■ P and y' — k ■ p. We refrain from stating explicitely 
the lengthy expressions for the quantities (<?o) y and (if' ) y here. 



2.2 Chebyshev Approximation 



The approximation of a function by Chebyshev polynomials of the first kind 
T n {z) is discussed in detail in ref. [22]. We briefly summarize the necessary 
formulae. Employing a convenient (albeit non-standard) normalization To = 
1 / y/2, the orthogonality relation reads, 

1 T n {z)T m {z) 7T 



-1 

Let {z k , k — 1 . . . K} denote the zeros of T K . The discrete orthogonality rela- 
tion is 

2 K 

~77 T m (z k )T n (z k ) = 5 mn . (42) 
A k=i 

We approximate a function F(p 2 ,p ■ P) by a finite sum of T's in the variable 
z = p ■ P: 



F(p 2 ,p-P)=Y,^ Fn (P 2 ) T n(z), (43) 

n=0 

^rriax H~l 

F"(p 2 ) = H)« J2 T n (z k )F(p 2 ,z k ). (44) 

k=l 



As before, the {z k } are the zeros of T„ max+1 . Note that for a finite n max the 
such defined Chebyshev moments F n (p 2 ) are not identical to the projections 
of eq. (43) using eq. (41), rather eqs. (43,44) are very close to the approx- 
imation of F by the minimax polynomial, which (among all polynomials of 
the same degree) has the smallest maximum deviation from the true function 
[22]. For z = z k and/or n max — > oo, projection and approximation are of course 
identical, courtesy of eq. (42). 

We use this method to approximate propagators and the exchange kernel in 
eqs. (39,40). Expanding the amplitudes of wave and vertex function as 
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Y z (p 2 ,z) = £ i m Yr(p 2 )T m (z), (45) 



m=0 



Y l {p\z)=Y J i n Y i n {p 2 )T n {z), (46) 



n=0 



and applying eqs. (43,44) to (g'o)^ in eq. (39), we obtain the matrix equation 



Y t n (p 2 ) = (9 r' nm (p 2 )Y J m (p 2 ), (47) 

O Pmax 

(^^(p 2 ) = — £ (^V, **) T B (z fc ) T m (z k ) . (48) 

Pmax ~r J- ^—i 

The projection onto the Yf was done using the orthogonality relation (41). 
In practice, we choose p max = n max + m max for a reliable approximation of 
the propagator matrix. The elements of the propagator matrix (goY^ nm are all 
real due to the explicit phase factor i n in the Chebyshev expansions (45,46). 

The Chebyshev approximation in the variables z and z' is now applied to the 
matrix (H') ij in eq. (40): 



+i ' l-m ax +i 

{H"f = {l-z' 2 ){H'f = £ E c^ st T s ^(z)T t ^(z') . (49) 

5 = 1 t=l 
2 2 m max + l n max + 1 

^"= (m +1)( „ +1) E E ^(*,«) 

x (H"r(p 2 ,k 2 ,y',z u ,z' v ). (50) 

The and {2^} are the zeros of the Chebyshev polynomial T minax+ i(z) and 
T nmax+ i(z'), respectively. We insert this expression into eq. (40), as well as the 
expansions (45,46) for wave and vertex function. After projecting onto the Y™ 
the integrations over z and z' on the r.h.s. are done using the orthogonality 
relation (41) and we obtain: 



+i ^max +i 

K"V) = 7 77 7 —^d\k\ dy' V V i n ~ m 

1 yP ' (m max +l)(n max + l)i 4tt 2 1 ' 7 ^ ^ 

x T m {z u )T n {z' v ) (H"y\p 2 ,k 2 ,y',z u ,z' v ) YJ>(k 2 ) . (51) 

Here the sum runs also over the label j and the Chebyshev moment label n. 

We have now succeeded to transform the original 4-dimensional integral equa- 
tion into a system of coupled one-dimensional equations. In summary, the 
system reads, 
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(52) 



Y t m (p 2 ) = J d\k\ H^' mn (k 2 ,p 2 ) Y?(k 2 ) , 



(53) 



o 



where indices appearing twice are summed over. Furthermore, the definition 



2.3 Solving the coupled one- dimensional integral equations 

The integration variable in eq. (53) is the absolute value of the four momentum 
k. It will be discretized on a mesh with nijM points, with typically n\u = 
20, ... , 50, and the integration is performed as a Gaussian quadrature. In the 
actual version of the program, the integration domain is first mapped onto the 
interval [-1,1] with the help of the transformation x = (\k\ — l)/(\k\ + 1) and 
then a Gauss-Legendre integration is used. 

After this step, the problem has become equivalent to finding the largest eigen- 
value of a matrix equation. To this end we employ an iterative method: Some 
initial values for the vertex functions, i.e. for the moments Y™, are chosen, 
and eq. (52) is used to calculate the corresponding expression for the wave 
function, i.e. for the moments Y™. Eq. (53) provides the moments Y™ in the 
next iteration step. If Y®(p 2 ) deviates more than a given accuracy (provided 
in the input file) for some momentum on the grid the updated moment Y®(p 2 ) 
at the grid point with the largest deviation is used to update the coupling 
constant g s . This value and renormalized functions Y™ are then the starting 
point for the next iteration step. Requiring an accuracy of the order 10~ 9 
convergence is usually reached after 20 - 30 iteration steps. 

The correct nucleon mass corresponds to the eigenvalue A = 1/g 2 . To de- 
termine it, in a first step eigenvalues A(Mi) and A(M 2 ) are calculated with 
Mi = max[m q ,m sc ,m ax ] and M 2 = mm[m g /r),m sc /(l-r)),m ax /(l-r))]. These 
starting values are motivated by the observation that for physically reasonable 
parameters the nucleon mass should be larger than the quark and diquark 
masses, and, for a given choice of 77, M 2 represents the upper limit beyond 
which singularities are encountered in the quark and diquark propagators. 
Now the nucleon mass M for which A(M) = 1/g 2 holds is determined by lin- 
ear bisection in the next step and by quadratic interpolation in all following 



H^ mn (k 2 ,p 2 ) 



(m max + l)(n 

max 

+ 1) 4tt 2 

x T m (z u )T n (z' v ) {H"f{p 2 



1 k 3 




■n—m 



(54) 



is used. 
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oet 


V 


m q 


m sc 




CO 


9s 


9a 


1V1 






[GeV] 


[GeV] 


[GeV] 


[GeV] 






[GeV] 


I 


0.36 


0.360 


0.6255 


0.6840 


0.95 


9.29 


6.97 


0.939 


II 


0.40 


0.425 


0.5977 


0.8314 


0.53 


22.10 


6.37 


0.939 



Table 1 

The two parameter sets together with the values of couplings and the bound mass 
M that arise for these sets. The maximal number of Chebyshev polynomials for 
vertex and wave function is given by m max = 8 and n max = 12 and the number 
of momentum grid points is n\ k \ = 40. The y' integration was done using Gaussian 
quadrature with n y = 32 grid points. 

steps. Typically 5-8 steps are needed to arrive at an accuracy of the order 
of 10~ 5 . During each step the integral equation has to be solved anew for 
its largest eigenvalue, thus illustrating the need for a fast algorithm which is 
provided here. 

As the system under consideration is a linear integral equation, the normal- 
ization for the functions Y™ and stays undetermined. Without loss of 
generality we fix the value of the lowest Chebychev moment of the first up- 
per component of the vertex function at the smallest absolute value on the 
momentum grid, 

n°(Pmin « 0) = 1, (55) 

and renormalize accordingly all functions Y™ and Y™ in the output routine. 
These vertex and wave functions are then the final result of the algorithm. 



3 Numerical Results 

To keep this article reasonably self-contained we will present numerical results 
only for two example sets of parameters. Quite a number of results obtained 
with this program can be found in refs. [9,13-15]. Altogether there are four 
parameters: the quark mass m q , the diquark masses m sc and m ax , and the 
width c of the dipole-shaped diquark-quark vertex (14). Diquark masses and 
the width determine the coupling constants g a and g s . Note that in actual 
applications g a and therefore m ax is determined by fitting the mass of the A 
baryon, see e.g. ref. [9] for more details. The equation for the A baryon is 
solved by an algorithm completely analogous to the one presented here. The 
results for two characteristic parameter sets are given in table 1. 
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Nucleon s waves Nucleon s waves 

Set I, zeroth Chebyshev moment Set II, zeroth Chebyshev moment 




0.05 0.1 0.15 0.2 0.25 0.3 0.05 0.1 0.15 0.2 0.25 0.3 



Ipl [GeV] |p| [GeV] 

Fig. 2. The leading Chebyshev moments of the functions Si, A2 and 
(1/3)^3 + (2/3) A5 related to the nucleon s waves. All functions are normalized 
by the condition S® = 1. 

Although the wave and the vertex function are no physical observables they 
do enter observable matrix elements (see e.g. refs. [13-15]) and therefore the 
strengths of the single components give a hint on their effect on observables. 
We have plotted the leading Chebyshev moments of the scalar functions de- 
scribing the nucleon s waves in figure 2. These are S® and (1/3) A3 + (2/3) A", 
describing the s waves for scalar and axialvector diquark and A\ that is con- 
nected with the virtual time component of the latter. We see for both sets that 
the functions A° are suppressed by a factor of 10 3 compared to the dominating 
5°. The strength of the other s wave associated with the axialvector diquark 
is roughly proportional to the ratio g a j 1 g s for the respective parameter set. 

In figure 3 we have plotted the th , 2 nd and 4 th Chebychev moments (the 
odd ones have been left out for clarity of the presentation) of the linear com- 
binations of the functions Si and M which belong to the eigenfunctions to 
the (3 quark) spin operator and the operator of angular momentum between 
quark and diquark [15]. One sees easily that in the nucleon rest frame higher 
Chebychev moments become unimportant very fast. 



4 Summary and Applications 

In summary, we have succeeded in developing an efficient algorithm for the 
calculation of the nucleon amplitudes in the seperable approximation to the 
Poincare covariant three-quark Faddeev equation. To this end we employ an 
approximate 0(4) symmetry of the resulting effective di quark-quark Bethe- 
Salpeter equation by using a Chebychev expansion. This allows to transform 
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Nucleon: 1=0, s=1/2 - Scalar Diquark 

S, 




0.1 0.2 0.3 

IPI [GeV] 

Nucleon: l=0, s=1/2 - Axialvector Diquark 

K 



- th Cheb. moment 
2 nd Cheb. moment 
4 ,h Cheb. moment 



0.1 0.2 



0.3 0.4 



IPl [GeV] 

Nucleon: l=0, s=1/2 - Axialvector Diquark 

1 / 3 A 3 + 2 / 3 a 5 




Cheb. moment 
2" Cheb. moment - 
4 th Cheb. moment 



0.1 0.2 0.3 0.4 0.5 

IPl [GeV] 

Nucleon: 1=2, s=3/2 - Axialvector Diquark 









,h Cheb. moment 




2"" Cheb. moment 




4 ,h Cheb. moment 



Nucleon: 1=1, s=1/2 -Scalar Diquark 







0* 


Cheb. moment - 






2 r 


Cheb. moment 






4' 


Cheb. moment 
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/ 

I / 
i / 









0.1 0.2 0.3 0.4 0.5 

IPl [GeV] 

Nucleon: 1=1 , s=1/2 - Axialvector Diquark 

K 









\\~~ 








o tB 


Cheb. moment 




2" 


Cheb. moment 




4* 


Cheb. moment 



0.1 0.2 0.3 0.4 0.5 

IPl [GeV] 

Nucleon: 1=1, s=1/2 - Axialvector Diquark 





th Cheb. moment 




2 rd Cheb. moment 




4 th Cheb. moment 











IPl [GeV] 

Nucleon: 1=1 , s=3/2 - Axialvector Diquark 

K-K 



0.2 0.3 

IPl [GeV] 




IPl [GeV] 



Fig. 3. Chebyshev moments of the scalar functions describing the strengths of 
quark-diquark partial waves with quantum number s for the total quark-diquark 
spin and quantum number I for the orbital angular momentum. 
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the full 4-dimensional problem into a manageable form. 

The chosen expansion in Dirac space into partial waves has been done covari- 
antly and the results for the partial wave strengths identify the leading and 
subleading components. The Chebyshev expansion for the scalar functions 
which describe the partial wave strengths has proved to be quite efficient and 
accurate, furthermore the expansion has been done for a Lorentz invariant 
variable, with the consequence that the boost of the wave and vertex function 
solutions is rendered feasible. 

The latter point is of utter importance because the applications of the pre- 
sented results to calculate physical observables require to Lorentz boost the 
Bethe-Salpeter wave function. Please note also that in Lorentz frames where 
the nucleon is moving higher Chebychev moments will become increasingly 
important with increasing nucleon momentum. Nevertheless, based on the 
presented algorithm nucleon form factors have been calculated up to several 
GeV of transfered momentum [14], and even production processes have been 
studied for a similar energy range [13]. 

Given the fact that the required CPU time to run this program is very mod- 
erate even on PCs, and that the underlying Poincare invariant nucleon model 
has a much broader application range than non-relativistic models we are con- 
fident that the presented model and the related algorithm will find wide-spread 
use in future investigations of nucleon properties. 



5 Description of the Program 

5. 1 The main program 

Before the main program we have included a module global to provide a 
global definition of the input parameters. Furthermore, the integration points 
and weights are defined globally. 

The main program calls the subroutine read_parameters, allocates mem- 
ory space for the used fields, and calls then the subroutine int_weights three 
times determining thereby the integration points and weights for the one angle 
and the one momentum integration. Note that the momenta pp G [0, oo] are 
mapped to the interval [0,1] for a Gauss-Legendre quadrature of the momen- 
tum integration. 

The subroutine dqnorms is called to calculate g s and g a as prescribed by 
eqs. (16,17). 



20 



Then the iteration loop for the nucleon mass is started with two initial values 
for the bound state mass M. During each iteration step the main algorithm, 
given in the subroutine bsemain is called which calculates the vertex function, 
the wave function and the coupling g s = l/y/X which corresponds to the mass 
M. The mass is updated by bisection until satisfactory agreement between 
1 / V\ and the scalar diquark normalization is reached. 

Calling the subroutine output concludes the main program. 

5.2 Input and Output 

Reading the input and writing out the vertex function, the wave function and 
the coupling g s is defered to the subroutines reacLparameters and output. 
The format in which the parameters are read in and the results are written 
out is obvious from the appended sample input and sample output. 

5.3 The subroutines "bsemain" and "bse_kernel" 

The first action to be taken in the subroutine bsemain is to call the subrou- 
tine bse_kernel to generate the kernel of the equation as given in eq. (54). 
This includes the integration over y' and the summation over the zeros of the 
Chebychev polynomials. A complex auxiliary array s is introduced whose real 
part, denoted r, is then used to provide the explicit expression for the kernel. 

Next, the propagator matrix (goY^ mn is generated, via the intermediary steps 
of specifying the numerator of the matrix (g'o)^ and the application of eq. (48). 

Then the Chebychev moments Y™ are initialized, all by the same function 
pe~ p . The iteration is implemented by a DO WHILE loop. In the iteration 
loop, first the Chebychev expansion of the propagators (denominators and 
numerators) and in a subsequent step the one of the product of quark and 
diquark propagators is determined. This allows to calculate the functions Y™, 
i.e. the wave function \l>, from eq. (20). Finally, eq. (19) with the calculated 
kernel is used to update the vertex function $ (the Y t m ). The new function 
F]° is used as described in sect. 3 to update g s and to test for convergence. 

5.4 Further Subroutines and Functions 

The subroutine int_weights calculates the integration points and weights for 
Gauss-Legendre quadrature. 
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A 

4 





o 


1 n 
1U 







0.8910 


U.ylob 


n no 1 o 


U.yzzo 


U.92ZO 


u.9zzy 


2 


91 93 




9352 


9360 


9362 


9363 


4 




0.9365 


0.9381 


0.9386 


0.9387 


0.9388 


6 






0.9387 


0.9389 


0.9390 


0.9390 


8 








0.9390 


0.9390 


0.9390 



Table 2 



Convergence of the bound state mass M [GeV] in terms of m max and n max . Due 
to weaker convergence of the wave function, only n max > m max is considered. For 
momentum and angular integration = n y = 20 grid points have been used. 



The RECURSIVE FUNCTION Chebyone determines the values of the Cheby- 
chev polynomials of the first kind. 

The FUNCTION dq implements the diquark function (14). 



6 Testing the program 



Of course, trivial tests establishing the independence of the number of inte- 
gration points etc. have been performed. We also verified that the vertex and 
wave functions are independent of the initializing functions. 

In table 2 we demonstrate the convergence of the bound state mass in terms of 
m max and n max . We have chosen parameter set I from table 1 for that purpose. 
Since in the scalar diquark channel the binding energy is just 46 MeV, one 
might not expect good convergence due to the nearby poles in quark and 
diquark propagators. Nevertheless, already for m max = 6 and n max = 8 the 
nucleon mass is exact up to 0.1 MeV. 

A very non-trivial test consists in changing the Mandelstam parameter 77, see 
the discussion the end of sect. 1. Please note that changing rj may necessi- 
tate an increase in the order of the Chebychev expansion in order to obtain 
comparable accuracy. 
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TEST RUN 



input file: 



quark mass : 0.36 

scalar diquark mass : 0.6255 

axialvector diquark mass : 0.684 

diquark width : 0.95 

type of diquark amplitude : 1 

exponent of diquark amplitude : 2 

momentum partitioning eta : 0.36 

npe - momentum grid |k|,|p| : 20 

ny - angle grid \vec k \cdot \vec p : 20 

mmax - Chebyshev accuracy VF : 6 

nmax - Chebyshev accuracy WF : 8 

accuracy for eigenvalue : l.D-9 

accuracy for mass : l.D-5 

output file : sa.out 



standard output: 



gs : 9.289520 ga : 6.965978 

mass iteration # 1 

M_N = 0.684000 
calculating H_kernel... 
generated VF Chebyshev moments 

1 2 3 4 5 6 H_kernel done! 

# it 1/sqrt (lambda) 

5 0.15781675145113D+02 
10 0.16128853436214D+02 
15 0.16125685785706D+02 
20 0.16125712415957D+02 
25 0.16125712185873D+02 
25 0.16125712185873D+02 
I gs-l/sqrt (lambda) I / gs: 0.735904 

mass iteration # 2 

M_N = 0.976984 
calculating H_kernel... 
generated VF Chebyshev moments 

1 2 3 4 5 6 H_kernel done! 

# it 1/sqrt (lambda) 

5 . 66709436074472D+01 
10 0.71005125223069D+01 
15 0.70958216340534D+01 
20 0.70958710453183D+01 
25 0.70958705246610D+01 
27 0.70958705292115D+01 
I gs-l/sqrt (lambda) I / gs: 0.236142 

mass iteration # 3 

M_N = 0.905808 
calculating H_kernel... 
generated VF Chebyshev moments 

1 2 3 4 5 6 H_kernel done! 

# it 1/sqrt (lambda) 
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5 0.99728969409881D+01 
10 0.10508348360659D+02 
15 0.10500476420175D+02 
20 0.10500595152654D+02 
25 0.10500593361326D+02 
29 0.10500593387021D+02 
I gs-l/sqrt (lambda) I / gs: 0.130370 

mass Iteration # 4 

M_N = 0.936576 
calculating H_kernel... 
generated VF Chebyshev moments 

1 2 3 4 5 6 H_kernel done! 

# it 1/sqrt (lambda) 

5 0.88748406217450D+01 
10 0.93909382462111D+01 
15 0.93830314492941D+01 
20 0.93831520978449D+01 
25 0.93831502564257D+01 
29 0.93831502831329D+01 
I gs-l/sqrt (lambda) I / gs: 0.010079 

mass iteration # 5 

M_N = 0.938826 
calculating H_kernel... 
generated VF Chebyshev moments 

1 2 3 4 5 6 H_kernel done! 

# it 1/sqrt (lambda) 

5 0.87863777879217D+01 
10 0.93005927655577D+01 
15 0.92927115567109D+01 
20 0.92928316413400D+01 
25 0.92928298111699D+01 
29 0.92928298376773D+01 
I gs-l/sqrt (lambda) I / gs: 0.000356 

mass iteration # 6 

M_N = 0.938908 
calculating H_kernel... 
generated VF Chebyshev moments 

1 2 3 4 5 6 H_kernel done! 

# it 1/sqrt (lambda) 

5 0.87831354990374D+01 
10 0.92972802606202D+01 
15 0.92894000862408D+01 
20 0.92895201477319D+01 
25 0.92895183180267D+01 
29 0.92895183445258D+01 
I gs-l/sqrt (lambda) I / gs: 0.000000 

done . . . 



output file: 

################################################################## 

# mq : 0.360000 ms : 0.625500 ma : 0.684000 

# M : 0.938908 ga/gs : 0.749875 gs : 9.289518 

# eta : 0.360000 
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# Dq . ampl . 

# mmax 

# npe 



RATIONAL co 

6 nmax 
20 ny 



0.950000 nexp 



20 



################################################################## 

# VERT 

# S_l 



# S_2 



X 


FUNCTION 
Ipl 









1 






2 






3 






4 






5 






6 





3448D-02 





1000D+01 





1038D-01 





7696D- 


-04 





6986D- 


-06 





7131D- 


-08 





8341D- 


-10 





1139D- 


-11 





1834D-01 





9948D+00 





5474D-01 





2145D- 


-02 





1029D- 


-03 





5536D- 


-05 





3407D- 


-06 





2456D- 


-07 





4590D-01 





9678D+00 





1305D+00 





1237D- 


-01 





1432D- 
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